source("scripts/recursos.R")
source("scripts/funciones.R")5 Diagnósticos de Residuos
5.1 Pretratamientos de la Información
Cargar Recursos Externos
Definición de variables generales
## geographic projections
crs_latlon <- "+proj=longlat +datum=WGS84 +no_defs"
crs_utm <-
"+proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"Definición de variables específicas
# List of crimes/delitos and cities/ciudades
crimes = c(
"Minor property offenses",
"Robberies",
"Burglaries",
"Drunkennes, damages and disorders",
"Family violence and aggression",
"Injuries, drugs and weapons",
"High violence and murder"
)
delitos = c("hurto", "robo", "asalto", "amen", "abus", "crim", "viol")
#ciudad de ejemplo
cities <- c("Coquimbo-Serena")
ciudades <- c("urb_4_18")
obs <- 1
delito <- delitos[1]Lectura de Datos Espaciales
Para el caso del presente ejemplo práctico se utiizará la ciudad de Coquimbo-Serena.
city <- ciudades[obs]
cityp <- readRDS(paste0("data/",city,".rds"))
cityp <- spTransform(cityp, CRS(crs_utm))El objeto espacial cityp cuyos puntos corresponde donde ha ucurrido eventos deictuales, cuyo formato SpatialPointDataFrame.
Matriz de vecindad espacial: Se utilizó la función propia spatial_weights() @(sec-spatial-weights)
# matriz de vecindad espacial
nb <- spatial_weights(city_df = cityp, nvec= 12)Transformaciones a los datos.
Se utilizó funciónes de transformacion general mencionadas en Section 6.2 que corresponden a las siguientes:
- Transformación Cúbica
- Transformación Logarítmica
- Cáculo del Lag Espacial
cityp <- cityp %>%
st_as_sf(crs = 4326) %>%
mutate(across(.cols = all_of(delitos),
.fns = t_cub, .names = "cub{.col}")) %>%
mutate(across(.cols = all_of(delitos),
.fns = t_log1p, .names = "log{.col}")) %>%
mutate(across(.cols = all_of(delitos),
.fns = lag_spatial, nb, .names = "lag{.col}")) %>%
as("Spatial")| zona | id | x | y | pflot | viol | crim | abus | amen | robo | hurto | asalto | denspob | educ | tamhog | constr | ofcom | vulzon | mobil | lpflot | neduc | ldenspob | lofcom | lconstr | cubhurto | cubrobo | cubasalto | cubamen | cubabus | cubcrim | cubviol | loghurto | logrobo | logasalto | logamen | logabus | logcrim | logviol | laghurto | lagrobo | lagasalto | lagamen | lagabus | lagcrim | lagviol |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 4101161005 | 22 | -71.239 | -29.862 | 1.000 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 144.756 | 10.867 | 3.170 | 4328.618 | 0.000 | 0.367 | 0.079 | 0.178 | 0.477 | 0.767 | 0.000 | 0.840 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0.000 | 0.000 | 0.000 | 0 | 0 | 0.083 | 0.083 | 0.083 | 0.417 | 0.250 | 0.083 | 0 |
| 4101161005 | 32 | -71.243 | -29.863 | 1.000 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 227.772 | 10.716 | 3.393 | 3830.009 | 0.000 | 0.367 | 0.079 | 0.178 | 0.466 | 0.846 | 0.000 | 0.829 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0.000 | 0.000 | 0.000 | 0 | 0 | 0.000 | 0.083 | 0.083 | 0.917 | 0.583 | 0.083 | 0 |
| 4101161005 | 35 | -71.240 | -29.863 | 1.034 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 171.133 | 11.801 | 3.172 | 4328.618 | 0.000 | 0.367 | 0.079 | 0.182 | 0.550 | 0.796 | 0.000 | 0.840 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0.693 | 0.000 | 0.000 | 0 | 0 | 0.000 | 0.083 | 0.000 | 0.500 | 0.250 | 0.083 | 0 |
| 4101161005 | 36 | -71.239 | -29.863 | 1.000 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 154.603 | 10.870 | 3.135 | 4328.618 | 0.000 | 0.367 | 0.079 | 0.178 | 0.478 | 0.779 | 0.000 | 0.840 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0.000 | 0.000 | 0.000 | 0 | 0 | 0.083 | 0.083 | 0.083 | 0.417 | 0.250 | 0.083 | 0 |
| 4101161005 | 37 | -71.238 | -29.863 | 1.025 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 155.145 | 10.858 | 3.148 | 3230.309 | 0.000 | 0.367 | 0.079 | 0.181 | 0.477 | 0.779 | 0.000 | 0.815 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0.000 | 0.000 | 0.000 | 0 | 0 | 0.083 | 0.083 | 0.083 | 0.000 | 0.250 | 0.083 | 0 |
| 4101161001 | 48 | -71.244 | -29.864 | 1.125 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 168.179 | 11.958 | 3.336 | 2384.222 | 14.682 | 0.292 | 0.122 | 0.193 | 0.562 | 0.793 | 0.574 | 0.788 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0.000 | 0.000 | 0.000 | 0 | 0 | 0.000 | 0.000 | 0.250 | 0.500 | 0.417 | 0.083 | 0 |
| 4101161005 | 49 | -71.243 | -29.864 | 1.071 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 212.743 | 10.978 | 3.415 | 3795.759 | 0.000 | 0.367 | 0.079 | 0.187 | 0.486 | 0.834 | 0.000 | 0.828 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0.000 | 0.000 | 0.000 | 0 | 0 | 0.000 | 0.083 | 0.250 | 0.917 | 0.583 | 0.083 | 0 |
| 4101161005 | 50 | -71.242 | -29.864 | 1.000 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 220.188 | 10.705 | 3.409 | 4026.613 | 0.000 | 0.367 | 0.079 | 0.178 | 0.465 | 0.840 | 0.000 | 0.834 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0.000 | 0.693 | 0.693 | 0 | 0 | 0.000 | 0.083 | 0.000 | 1.000 | 0.417 | 0.083 | 0 |
| 4101161005 | 51 | -71.241 | -29.864 | 1.021 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 243.501 | 10.120 | 3.337 | 4328.618 | 0.000 | 0.367 | 0.079 | 0.180 | 0.419 | 0.857 | 0.000 | 0.840 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0.000 | 0.693 | 0.693 | 0 | 0 | 0.000 | 0.083 | 0.083 | 0.667 | 0.250 | 0.000 | 0 |
| 4101161005 | 52 | -71.240 | -29.864 | 1.022 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 320.023 | 10.720 | 3.324 | 4328.618 | 0.000 | 0.367 | 0.079 | 0.180 | 0.466 | 0.905 | 0.000 | 0.840 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0.000 | 0.000 | 0.000 | 0 | 0 | 0.000 | 0.083 | 0.167 | 0.583 | 0.250 | 0.000 | 0 |
5.2 Evaluación De Residuos
5.2.1 Modelos MLB
Lectura del Modelo.
Warning in spTransform(xSP, CRSobj, ...): NULL source CRS comment, falling back
to PROJ string
Warning in wkt(obj): CRS object has no comment
model_mlb <- readRDS(paste0("output/mlb_",city,"_",delito,".rds"))5.2.1.1 Extracción de Residuos
cityp$res_fit_mlb <- residuals(model_mlb)
cityp$fitted_fit_mlb <- fitted(model_mlb)
cityp$sd_breaks_mlb <- scale(cityp$res_fit_mlb)[,1]breaks_qt <- classInt::classIntervals(cityp$sd_breaks_mlb,
n = 5, style = "fisher")
my_breaks <- round(breaks_qt$brks,2)city_sf <- cityp %>% st_as_sf(crs = 32719)
residuos_mv <- mapview(city_sf,zcol="res_fit_mlb", at =my_breaks, cex=2)
residuos_mv5.2.1.2 Global Moran
Global Moran a Variable Depediente
spdep::moran.test(x = cityp$laghurto, listw = nb)
Moran I test under randomisation
data: cityp$laghurto
weights: nb
Moran I statistic standard deviate = 174.4, p-value <
0.00000000000000022
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.7640052224 -0.0001215805 0.0000191963
# moran.plot(x =cityp$laghurto, listw = nb, labels=as.character(cityp$id))Global Moran de Residuos
spdep::moran.test(x = cityp$res_fit_mlb, listw = nb)
Moran I test under randomisation
data: cityp$res_fit_mlb
weights: nb
Moran I statistic standard deviate = 19.256, p-value <
0.00000000000000022
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.08478830236 -0.00012158055 0.00001944476
# moran.plot(x =cityp$res_fit_mlb, listw = nb, labels=as.character(cityp$id))5.2.1.3 Local Moran de Residuos
# Calcular Local Moran
lmoran = localmoran(cityp$res_fit_mlb, nb)
cityp$res_scaled = as.numeric(scale(cityp$res_fit_mlb))
cityp$lag_scaled = lag.listw(nb, cityp$res_scaled)
cityp@data = cbind(cityp@data, lmoran = as.data.frame(lmoran)[, 5])
# Umbral de significancia estadistica
pval=0.05
# Definir cuadrantes
cityp[(cityp$res_scaled >= 0 & cityp$lag_scaled >= 0) & (cityp$lmoran <= pval), "clusterM"] = "HH" # plot
cityp[(cityp$res_scaled <= 0 & cityp$lag_scaled <= 0) & (cityp$lmoran <= pval), "clusterM"] = "LL" # plot
cityp[(cityp$res_scaled >= 0 & cityp$lag_scaled <= 0) & (cityp$lmoran <= pval), "clusterM"] = "HL"
cityp[(cityp$res_scaled <= 0 & cityp$lag_scaled >= 0) & (cityp$lmoran <= pval), "clusterM"] = "LH"
cityp[(lmoran[, 5] > pval), "clusterM"] = "NS"
table(cityp$clusterM)
HH HL LH LL NS
284 33 295 331 7283
city_sf <- cityp %>% st_as_sf(crs = 32719) %>%
filter(clusterM != "NS")
residuos_moran <- mapview(city_sf,zcol="clusterM", cex=2 )
residuos_moran5.2.2 Modelos SAR
Warning in spTransform(xSP, CRSobj, ...): NULL source CRS comment, falling back
to PROJ string
Warning in wkt(obj): CRS object has no comment
model_sar <- readRDS(paste0("output/sar_",city,"_",delito,".rds"))5.2.2.1 Extracción de Residuos SAR
cityp$res_fit_sar <- residuals(model_sar)
cityp$fitted_fit_sar <- fitted(model_sar)This method assumes the response is known - see manual page
cityp$sd_breaks_sar <- scale(cityp$res_fit_sar)[,1]breaks_qt <- classInt::classIntervals(cityp$sd_breaks_sar,
n = 5, style = "fisher")
my_breaks <- round(breaks_qt$brks,2)city_sf <- cityp %>% st_as_sf(crs = 32719)
residuos_mv <- mapview(city_sf,zcol="res_fit_sar", at =my_breaks, cex=2)
residuos_mv5.2.2.2 Global Moran
Global Moran a Variable Depediente
spdep::moran.test(x = cityp$cubhurto, listw = nb)
Moran I test under randomisation
data: cityp$cubhurto
weights: nb
Moran I statistic standard deviate = 79.737, p-value <
0.00000000000000022
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.35153153375 -0.00012158055 0.00001944977
# moran.plot(x =cityp$laghurto, listw = nb, labels=as.character(cityp$id))Global Moran de Residuos
spdep::moran.test(x = cityp$res_fit_sar, listw = nb)
Moran I test under randomisation
data: cityp$res_fit_sar
weights: nb
Moran I statistic standard deviate = -1.2062, p-value = 0.8861
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
-0.00544188551 -0.00012158055 0.00001945508
# moran.plot(x =cityp$res_fit_mlb, listw = nb, labels=as.character(cityp$id))5.2.2.3 Local Moran de Residuos
# Calcular Local Moran
lmoran = localmoran(cityp$res_fit_sar, nb)
cityp$res_scaled = as.numeric(scale(cityp$res_fit_sar))
cityp$lag_scaled = lag.listw(nb, cityp$res_scaled)
cityp@data = cbind(cityp@data, lmoran = as.data.frame(lmoran)[, 5])
# Umbral de significancia estadistica
pval=0.05
# Definir cuadrantes
cityp[(cityp$res_scaled >= 0 & cityp$lag_scaled >= 0) & (cityp$lmoran <= pval), "clusterM"] = "HH" # plot
cityp[(cityp$res_scaled <= 0 & cityp$lag_scaled <= 0) & (cityp$lmoran <= pval), "clusterM"] = "LL" # plot
cityp[(cityp$res_scaled >= 0 & cityp$lag_scaled <= 0) & (cityp$lmoran <= pval), "clusterM"] = "HL"
cityp[(cityp$res_scaled <= 0 & cityp$lag_scaled >= 0) & (cityp$lmoran <= pval), "clusterM"] = "LH"
cityp[(lmoran[, 5] > pval), "clusterM"] = "NS"
table(cityp$clusterM)
HH HL LH LL NS
214 19 156 32 7805
city_sf <- cityp %>% st_as_sf(crs = 32719) %>%
filter(clusterM != "NS")
residuos_moran <- mapview(city_sf,zcol="clusterM", cex=2 )
residuos_moran